home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / implicit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  24.0 KB  |  815 lines
  1. /*
  2.  * C code from the article
  3.  * "An Implicit Surface Polygonizer"
  4.  * by Jules Bloomenthal, jbloom@beauty.gmu.edu
  5.  * in "Graphics Gems IV", Academic Press, 1994
  6.  */
  7.  
  8. /* implicit.c
  9.  *     an implicit surface polygonizer, translated from Mesa
  10.  *     applications should call polygonize()
  11.  *
  12.  * to compile a test program for ASCII output:
  13.  *     cc implicit.c -o implicit -lm
  14.  *
  15.  * to compile a test program for display on an SGI workstation:
  16.  *     cc -DSGIGFX implicit.c -o implicit -lgl_s -lm
  17.  *
  18.  * Authored by Jules Bloomenthal, Xerox PARC.
  19.  * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
  20.  * Permission is granted to reproduce, use and distribute this code for
  21.  * any and all purposes, provided that this notice appears in all copies. */
  22.  
  23. #include <math.h>
  24. #include <stdio.h>
  25. #include <sys/types.h>
  26.  
  27. #define TET    0  /* use tetrahedral decomposition */
  28. #define NOTET    1  /* no tetrahedral decomposition  */
  29.  
  30. #define RES    10 /* # converge iterations    */
  31.  
  32. #define L    0  /* left direction:    -x, -i */
  33. #define R    1  /* right direction:    +x, +i */
  34. #define B    2  /* bottom direction: -y, -j */
  35. #define T    3  /* top direction:    +y, +j */
  36. #define N    4  /* near direction:    -z, -k */
  37. #define F    5  /* far direction:    +z, +k */
  38. #define LBN    0  /* left bottom near corner  */
  39. #define LBF    1  /* left bottom far corner   */
  40. #define LTN    2  /* left top near corner     */
  41. #define LTF    3  /* left top far corner      */
  42. #define RBN    4  /* right bottom near corner */
  43. #define RBF    5  /* right bottom far corner  */
  44. #define RTN    6  /* right top near corner    */
  45. #define RTF    7  /* right top far corner     */
  46.  
  47. /* the LBN corner of cube (i, j, k), corresponds with location
  48.  * (start.x+(i-.5)*size, start.y+(j-.5)*size, start.z+(k-.5)*size) */
  49.  
  50. #define RAND()        ((rand()&32767)/32767.)    /* random number between 0 and 1 */
  51. #define HASHBIT        (5)
  52. #define HASHSIZE    (size_t)(1<<(3*HASHBIT))   /* hash table size (32768) */
  53. #define MASK        ((1<<HASHBIT)-1)
  54. #define HASH(i,j,k) ((((((i)&MASK)<<HASHBIT)|((j)&MASK))<<HASHBIT)|((k)&MASK))
  55. #define BIT(i, bit) (((i)>>(bit))&1)
  56. #define FLIP(i,bit) ((i)^1<<(bit)) /* flip the given bit of i */
  57.  
  58. typedef struct point {           /* a three-dimensional point */
  59.     double x, y, z;           /* its coordinates */
  60. } POINT;
  61.  
  62. typedef struct test {           /* test the function for a signed value */
  63.     POINT p;               /* location of test */
  64.     double value;           /* function value at p */
  65.     int ok;               /* if value is of correct sign */
  66. } TEST;
  67.  
  68. typedef struct vertex {           /* surface vertex */
  69.     POINT position, normal;       /* position and surface normal */
  70. } VERTEX;
  71.  
  72. typedef struct vertices {       /* list of vertices in polygonization */
  73.     int count, max;           /* # vertices, max # allowed */
  74.     VERTEX *ptr;           /* dynamically allocated */
  75. } VERTICES;
  76.  
  77. typedef struct corner {           /* corner of a cube */
  78.     int i, j, k;           /* (i, j, k) is index within lattice */
  79.     double x, y, z, value;       /* location and function value */
  80. } CORNER;
  81.  
  82. typedef struct cube {           /* partitioning cell (cube) */
  83.     int i, j, k;           /* lattice location of cube */
  84.     CORNER *corners[8];           /* eight corners */
  85. } CUBE;
  86.  
  87. typedef struct cubes {           /* linked list of cubes acting as stack */
  88.     CUBE cube;               /* a single cube */
  89.     struct cubes *next;           /* remaining elements */
  90. } CUBES;
  91.  
  92. typedef struct centerlist {       /* list of cube locations */
  93.     int i, j, k;           /* cube location */
  94.     struct centerlist *next;       /* remaining elements */
  95. } CENTERLIST;
  96.  
  97. typedef struct cornerlist {       /* list of corners */
  98.     int i, j, k;           /* corner id */
  99.     double value;           /* corner value */
  100.     struct cornerlist *next;       /* remaining elements */
  101. } CORNERLIST;
  102.  
  103. typedef struct edgelist {       /* list of edges */
  104.     int i1, j1, k1, i2, j2, k2;       /* edge corner ids */
  105.     int vid;               /* vertex id */
  106.     struct edgelist *next;       /* remaining elements */
  107. } EDGELIST;
  108.  
  109. typedef struct intlist {       /* list of integers */
  110.     int i;               /* an integer */
  111.     struct intlist *next;       /* remaining elements */
  112. } INTLIST;
  113.  
  114. typedef struct intlists {       /* list of list of integers */
  115.     INTLIST *list;           /* a list of integers */
  116.     struct intlists *next;       /* remaining elements */
  117. } INTLISTS;
  118.  
  119. typedef struct process {       /* parameters, function, storage */
  120.     double (*function)();       /* implicit surface function */
  121.     int (*triproc)();           /* triangle output function */
  122.     double size, delta;           /* cube size, normal delta */
  123.     int bounds;               /* cube range within lattice */
  124.     POINT start;           /* start point on surface */
  125.     CUBES *cubes;           /* active cubes */
  126.     VERTICES vertices;           /* surface vertices */
  127.     CENTERLIST **centers;       /* cube center hash table */
  128.     CORNERLIST **corners;       /* corner value hash table */
  129.     EDGELIST **edges;           /* edge and vertex id hash table */
  130. } PROCESS;
  131.  
  132. void *calloc();
  133. char *mycalloc();
  134.  
  135.  
  136. /**** A Test Program ****/
  137.  
  138.  
  139. /* torus: a torus with major, minor radii = 0.5, 0.1, try size = .05 */
  140.  
  141. double torus (x, y, z)
  142. double x, y, z;
  143. {
  144.     double x2 = x*x, y2 = y*y, z2 = z*z;
  145.     double a = x2+y2+z2+(0.5*0.5)-(0.1*0.1);
  146.     return a*a-4.0*(0.5*0.5)*(y2+z2);
  147. }
  148.  
  149.  
  150. /* sphere: an inverse square function (always positive) */
  151.  
  152. double sphere (x, y, z)
  153. double x, y, z;
  154. {
  155.     double rsq = x*x+y*y+z*z;
  156.     return 1.0/(rsq < 0.00001? 0.00001 : rsq);
  157. }
  158.  
  159.  
  160. /* blob: a three-pole blend function, try size = .1 */
  161.  
  162. double blob (x, y, z)
  163. double x, y, z;
  164. {
  165.     return 4.0-sphere(x+1.0,y,z)-sphere(x,y+1.0,z)-sphere(x,y,z+1.0);
  166. }
  167.  
  168.  
  169. #ifdef SGIGFX /**************************************************************/
  170.  
  171. #include "gl.h"
  172.  
  173. /* triangle: called by polygonize() for each triangle; set SGI lines */
  174.  
  175. triangle (i1, i2, i3, vertices)
  176. int i1, i2, i3;
  177. VERTICES vertices;
  178. {
  179.     float v[3];
  180.     int i, ids[3];
  181.     ids[0] = i1;
  182.     ids[1] = i2;
  183.     ids[2] = i3;
  184.     bgnclosedline();
  185.     for (i = 0; i < 3; i++) {
  186.     POINT *p = &vertices.ptr[ids[i]].position;
  187.     v[0] = p->x; v[1] = p->y; v[2] = p->z;
  188.     v3f(v);
  189.     }
  190.     endclosedline();
  191.     return 1;
  192. }
  193.  
  194.  
  195. /* main: call polygonize() with torus function
  196.  * display lines on SGI */
  197.  
  198. main ()
  199. {
  200.     char *err, *polygonize();
  201.  
  202.     keepaspect(1, 1);
  203.     winopen("implicit");
  204.     doublebuffer();
  205.     gconfig();
  206.     perspective(450, 1.0/1.0, 0.1, 10.0);
  207.     color(7);
  208.     clear();
  209.     swapbuffers();
  210.     makeobj(1);
  211.     if ((err = polygonize(torus, .05, 20, 0.,0.,0., triangle, TET)) != NULL) {
  212.     fprintf(stderr, "%s\n", err);
  213.     exit(1);
  214.     }
  215.     closeobj();
  216.     translate(0.0, 0.0, -2.0);
  217.     pushmatrix();
  218.     while(1) { /* spin the object */
  219.     reshapeviewport();
  220.     color(7);
  221.     clear();
  222.     color(0);
  223.     callobj(1);
  224.     rot(0.8, 'x');
  225.     rot(0.3, 'y');
  226.     rot(0.1, 'z');
  227.     swapbuffers();
  228.  
  229.     }
  230. }
  231.  
  232. #else /***********************************************************************/
  233.  
  234. int gntris;         /* global needed by application */
  235. VERTICES gvertices;  /* global needed by application */
  236.  
  237.  
  238. /* triangle: called by polygonize() for each triangle; write to stdout */
  239.  
  240. triangle (i1, i2, i3, vertices)
  241. int i1, i2, i3;
  242. VERTICES vertices;
  243. {
  244.     gvertices = vertices;
  245.     gntris++;
  246.     fprintf(stdout, "%d %d %d\n", i1, i2, i3);
  247.     return 1;
  248. }
  249.  
  250.  
  251. /* main: call polygonize() with torus function
  252.  * write points-polygon formatted data to stdout */
  253.  
  254. main ()
  255.     {
  256.     int i;
  257.     char *err, *polygonize();
  258.     gntris = 0;
  259.     fprintf(stdout, "triangles\n\n");
  260.     if ((err = polygonize(torus, .05, 20, 0.,0.,0., triangle, TET)) != NULL) {
  261.     fprintf(stdout, "%s\n", err);
  262.     exit(1);
  263.     }
  264.     fprintf(stdout, "\n%d triangles, %d vertices\n", gntris, gvertices.count);
  265.     fprintf(stdout, "\nvertices\n\n");
  266.     for (i = 0; i < gvertices.count; i++) {
  267.     VERTEX v;
  268.     v = gvertices.ptr[i];
  269.     fprintf(stdout, "%f  %f     %f\t%f     %f  %f\n",
  270.         v.position.x, v.position.y,     v.position.z,
  271.         v.normal.x,      v.normal.y,     v.normal.z);
  272.     }
  273.     fprintf(stderr, "%d triangles, %d vertices\n", gntris, gvertices.count);
  274.     exit(0);
  275. }
  276.  
  277. #endif /**********************************************************************/
  278.  
  279.  
  280. /**** An Implicit Surface Polygonizer ****/
  281.  
  282.  
  283. /* polygonize: polygonize the implicit surface function
  284.  *   arguments are:
  285.  *     double function (x, y, z)
  286.  *         double x, y, z (an arbitrary 3D point)
  287.  *         the implicit surface function
  288.  *         return negative for inside, positive for outside
  289.  *     double size
  290.  *         width of the partitioning cube
  291.  *     int bounds
  292.  *         max. range of cubes (+/- on the three axes) from first cube
  293.  *     double x, y, z
  294.  *         coordinates of a starting point on or near the surface
  295.  *         may be defaulted to 0., 0., 0.
  296.  *     int triproc (i1, i2, i3, vertices)
  297.  *         int i1, i2, i3 (indices into the vertex array)
  298.  *         VERTICES vertices (the vertex array, indexed from 0)
  299.  *         called for each triangle
  300.  *         the triangle coordinates are (for i = i1, i2, i3):
  301.  *         vertices.ptr[i].position.x, .y, and .z
  302.  *         vertices are ccw when viewed from the out (positive) side
  303.  *         in a left-handed coordinate system
  304.  *         vertex normals point outwards
  305.  *         return 1 to continue, 0 to abort
  306.  *     int mode
  307.  *         TET: decompose cube and polygonize six tetrahedra
  308.  *         NOTET: polygonize cube directly
  309.  *   returns error or NULL
  310.  */
  311.  
  312. char *polygonize (function, size, bounds, x, y, z, triproc, mode)
  313. double (*function)(), size, x, y, z;
  314. int bounds, (*triproc)(), mode;
  315. {
  316.     PROCESS p;
  317.     int n, noabort;
  318.     CORNER *setcorner();
  319.     TEST in, out, find();
  320.  
  321.     p.function = function;
  322.     p.triproc = triproc;
  323.     p.size = size;
  324.     p.bounds = bounds;
  325.     p.delta = size/(double)(RES*RES);
  326.  
  327.     /* allocate hash tables and build cube polygon table: */
  328.     p.centers = (CENTERLIST **) mycalloc(HASHSIZE,sizeof(CENTERLIST *));
  329.     p.corners = (CORNERLIST **) mycalloc(HASHSIZE,sizeof(CORNERLIST *));
  330.     p.edges =    (EDGELIST   **) mycalloc(2*HASHSIZE,sizeof(EDGELIST *));
  331.     makecubetable();
  332.  
  333.     /* find point on surface, beginning search at (x, y, z): */
  334.     srand(1);
  335.     in = find(1, &p, x, y, z);
  336.     out = find(0, &p, x, y, z);
  337.     if (!in.ok || !out.ok) return "can't find starting point";
  338.     converge(&in.p, &out.p, in.value, p.function, &p.start);
  339.  
  340.     /* push initial cube on stack: */
  341.     p.cubes = (CUBES *) mycalloc(1, sizeof(CUBES)); /* list of 1 */
  342.     p.cubes->cube.i = p.cubes->cube.j = p.cubes->cube.k = 0;
  343.     p.cubes->next = NULL;
  344.  
  345.     /* set corners of initial cube: */
  346.     for (n = 0; n < 8; n++)
  347.     p.cubes->cube.corners[n] = setcorner(&p, BIT(n,2), BIT(n,1), BIT(n,0));
  348.  
  349.     p.vertices.count = p.vertices.max = 0; /* no vertices yet */
  350.     p.vertices.ptr = NULL;
  351.  
  352.     setcenter(p.centers, 0, 0, 0);
  353.  
  354.     while (p.cubes != NULL) { /* process active cubes till none left */
  355.     CUBE c;
  356.     CUBES *temp = p.cubes;
  357.     c = p.cubes->cube;
  358.  
  359.     noabort = mode == TET?
  360.            /* either decompose into tetrahedra and polygonize: */
  361.            dotet(&c, LBN, LTN, RBN, LBF, &p) &&
  362.            dotet(&c, RTN, LTN, LBF, RBN, &p) &&
  363.            dotet(&c, RTN, LTN, LTF, LBF, &p) &&
  364.            dotet(&c, RTN, RBN, LBF, RBF, &p) &&
  365.            dotet(&c, RTN, LBF, LTF, RBF, &p) &&
  366.            dotet(&c, RTN, LTF, RTF, RBF, &p)
  367.            :
  368.            /* or polygonize the cube directly: */
  369.            docube(&c, &p);
  370.     if (! noabort) return "aborted";
  371.  
  372.     /* pop current cube from stack */
  373.     p.cubes = p.cubes->next;
  374.     free((char *) temp);
  375.     /* test six face directions, maybe add to stack: */
  376.     testface(c.i-1, c.j, c.k, &c, L, LBN, LBF, LTN, LTF, &p);
  377.     testface(c.i+1, c.j, c.k, &c, R, RBN, RBF, RTN, RTF, &p);
  378.     testface(c.i, c.j-1, c.k, &c, B, LBN, LBF, RBN, RBF, &p);
  379.     testface(c.i, c.j+1, c.k, &c, T, LTN, LTF, RTN, RTF, &p);
  380.     testface(c.i, c.j, c.k-1, &c, N, LBN, LTN, RBN, RTN, &p);
  381.     testface(c.i, c.j, c.k+1, &c, F, LBF, LTF, RBF, RTF, &p);
  382.     }
  383.     return NULL;
  384. }
  385.  
  386.  
  387. /* testface: given cube at lattice (i, j, k), and four corners of face,
  388.  * if surface crosses face, compute other four corners of adjacent cube
  389.  * and add new cube to cube stack */
  390.  
  391. testface (i, j, k, old, face, c1, c2, c3, c4, p)
  392. CUBE *old;
  393. PROCESS *p;
  394. int i, j, k, face, c1, c2, c3, c4;
  395. {
  396.     CUBE new;
  397.     CUBES *oldcubes = p->cubes;
  398.     CORNER *setcorner();
  399.     static int facebit[6] = {2, 2, 1, 1, 0, 0};
  400.     int n, pos = old->corners[c1]->value > 0.0 ? 1 : 0, bit = facebit[face];
  401.  
  402.     /* test if no surface crossing, cube out of bounds, or already visited: */
  403.     if ((old->corners[c2]->value > 0) == pos &&
  404.     (old->corners[c3]->value > 0) == pos &&
  405.     (old->corners[c4]->value > 0) == pos) return;
  406.     if (abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;
  407.     if (setcenter(p->centers, i, j, k)) return;
  408.  
  409.     /* create new cube: */
  410.     new.i = i;
  411.     new.j = j;
  412.     new.k = k;
  413.     for (n = 0; n < 8; n++) new.corners[n] = NULL;
  414.     new.corners[FLIP(c1, bit)] = old->corners[c1];
  415.     new.corners[FLIP(c2, bit)] = old->corners[c2];
  416.     new.corners[FLIP(c3, bit)] = old->corners[c3];
  417.     new.corners[FLIP(c4, bit)] = old->corners[c4];
  418.     for (n = 0; n < 8; n++)
  419.     if (new.corners[n] == NULL)
  420.         new.corners[n] = setcorner(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0));
  421.  
  422.     /*add cube to top of stack: */
  423.     p->cubes = (CUBES *) mycalloc(1, sizeof(CUBES));
  424.     p->cubes->cube = new;
  425.     p->cubes->next = oldcubes;
  426. }
  427.  
  428.  
  429. /* setcorner: return corner with the given lattice location
  430.    set (and cache) its function value */
  431.  
  432. CORNER *setcorner (p, i, j, k)
  433. int i, j, k;
  434. PROCESS *p;
  435. {
  436.     /* for speed, do corner value caching here */
  437.     CORNER *c = (CORNER *) mycalloc(1, sizeof(CORNER));
  438.     int index = HASH(i, j, k);
  439.     CORNERLIST *l = p->corners[index];
  440.     c->i = i; c->x = p->start.x+((double)i-.5)*p->size;
  441.     c->j = j; c->y = p->start.y+((double)j-.5)*p->size;
  442.     c->k = k; c->z = p->start.z+((double)k-.5)*p->size;
  443.     for (; l != NULL; l = l->next)
  444.     if (l->i == i && l->j == j && l->k == k) {
  445.         c->value = l->value;
  446.         return c;
  447.         }
  448.     l = (CORNERLIST *) mycalloc(1, sizeof(CORNERLIST));
  449.     l->i = i; l->j = j; l->k = k;
  450.     l->value = c->value = p->function(c->x, c->y, c->z);
  451.     l->next = p->corners[index];
  452.     p->corners[index] = l;
  453.     return c;
  454. }
  455.  
  456.  
  457. /* find: search for point with value of given sign (0: neg, 1: pos) */
  458.  
  459. TEST find (sign, p, x, y, z)
  460. int sign;
  461. PROCESS *p;
  462. double x, y, z;
  463. {
  464.     int i;
  465.     TEST test;
  466.     double range = p->size;
  467.     test.ok = 1;
  468.     for (i = 0; i < 10000; i++) {
  469.     test.p.x = x+range*(RAND()-0.5);
  470.     test.p.y = y+range*(RAND()-0.5);
  471.     test.p.z = z+range*(RAND()-0.5);
  472.     test.value = p->function(test.p.x, test.p.y, test.p.z);
  473.     if (sign == (test.value > 0.0)) return test;
  474.     range = range*1.0005; /* slowly expand search outwards */
  475.     }
  476.     test.ok = 0;
  477.     return test;
  478. }
  479.  
  480.  
  481. /**** Tetrahedral Polygonization ****/
  482.  
  483.  
  484. /* dotet: triangulate the tetrahedron
  485.  * b, c, d should appear clockwise when viewed from a
  486.  * return 0 if client aborts, 1 otherwise */
  487.  
  488. int dotet (cube, c1, c2, c3, c4, p)
  489. CUBE *cube;
  490. int c1, c2, c3, c4;
  491. PROCESS *p;
  492. {
  493.     CORNER *a = cube->corners[c1];
  494.     CORNER *b = cube->corners[c2];
  495.     CORNER *c = cube->corners[c3];
  496.     CORNER *d = cube->corners[c4];
  497.     int index = 0, apos, bpos, cpos, dpos, e1, e2, e3, e4, e5, e6;
  498.     if (apos = (a->value > 0.0)) index += 8;
  499.     if (bpos = (b->value > 0.0)) index += 4;
  500.     if (cpos = (c->value > 0.0)) index += 2;
  501.     if (dpos = (d->value > 0.0)) index += 1;
  502.     /* index is now 4-bit number representing one of the 16 possible cases */
  503.     if (apos != bpos) e1 = vertid(a, b, p);
  504.     if (apos != cpos) e2 = vertid(a, c, p);
  505.     if (apos != dpos) e3 = vertid(a, d, p);
  506.     if (bpos != cpos) e4 = vertid(b, c, p);
  507.     if (bpos != dpos) e5 = vertid(b, d, p);
  508.     if (cpos != dpos) e6 = vertid(c, d, p);
  509.     /* 14 productive tetrahedral cases (0000 and 1111 do not yield polygons */
  510.     switch (index) {
  511.     case 1:     return p->triproc(e5, e6, e3, p->vertices);
  512.     case 2:     return p->triproc(e2, e6, e4, p->vertices);
  513.     case 3:     return p->triproc(e3, e5, e4, p->vertices) &&
  514.             p->triproc(e3, e4, e2, p->vertices);
  515.     case 4:     return p->triproc(e1, e4, e5, p->vertices);
  516.     case 5:     return p->triproc(e3, e1, e4, p->vertices) &&
  517.             p->triproc(e3, e4, e6, p->vertices);
  518.     case 6:     return p->triproc(e1, e2, e6, p->vertices) &&
  519.             p->triproc(e1, e6, e5, p->vertices);
  520.     case 7:     return p->triproc(e1, e2, e3, p->vertices);
  521.     case 8:     return p->triproc(e1, e3, e2, p->vertices);
  522.     case 9:     return p->triproc(e1, e5, e6, p->vertices) &&
  523.             p->triproc(e1, e6, e2, p->vertices);
  524.     case 10: return p->triproc(e1, e3, e6, p->vertices) &&
  525.             p->triproc(e1, e6, e4, p->vertices);
  526.     case 11: return p->triproc(e1, e5, e4, p->vertices);
  527.     case 12: return p->triproc(e3, e2, e4, p->vertices) &&
  528.             p->triproc(e3, e4, e5, p->vertices);
  529.     case 13: return p->triproc(e6, e2, e4, p->vertices);
  530.     case 14: return p->triproc(e5, e3, e6, p->vertices);
  531.     }
  532.     return 1;
  533. }
  534.  
  535.  
  536. /**** Cubical Polygonization (optional) ****/
  537.  
  538.  
  539. #define LB    0  /* left bottom edge    */
  540. #define LT    1  /* left top edge    */
  541. #define LN    2  /* left near edge    */
  542. #define LF    3  /* left far edge    */
  543. #define RB    4  /* right bottom edge */
  544. #define RT    5  /* right top edge    */
  545. #define RN    6  /* right near edge    */
  546. #define RF    7  /* right far edge    */
  547. #define BN    8  /* bottom near edge    */
  548. #define BF    9  /* bottom far edge    */
  549. #define TN    10 /* top near edge    */
  550. #define TF    11 /* top far edge    */
  551.  
  552. static INTLISTS *cubetable[256];
  553.  
  554. /*            edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
  555. static int corner1[12]       = {LBN,LTN,LBN,LBF,RBN,RTN,RBN,RBF,LBN,LBF,LTN,LTF};
  556. static int corner2[12]       = {LBF,LTF,LTN,LTF,RBF,RTF,RTN,RTF,RBN,RBF,RTN,RTF};
  557. static int leftface[12]       = {B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F};
  558.                  /* face on left when going corner1 to corner2 */
  559. static int rightface[12]   = {L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T};
  560.                  /* face on right when going corner1 to corner2 */
  561.  
  562.  
  563. /* docube: triangulate the cube directly, without decomposition */
  564.  
  565. int docube (cube, p)
  566. CUBE *cube;
  567. PROCESS *p;
  568. {
  569.     INTLISTS *polys;
  570.     int i, index = 0;
  571.     for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0) index += (1<<i);
  572.     for (polys = cubetable[index]; polys; polys = polys->next) {
  573.     INTLIST *edges;
  574.     int a = -1, b = -1, count = 0;
  575.     for (edges = polys->list; edges; edges = edges->next) {
  576.         CORNER *c1 = cube->corners[corner1[edges->i]];
  577.         CORNER *c2 = cube->corners[corner2[edges->i]];
  578.         int c = vertid(c1, c2, p);
  579.         if (++count > 2 && ! p->triproc(a, b, c, p->vertices)) return 0;
  580.         if (count < 3) a = b;
  581.         b = c;
  582.     }
  583.     }
  584.     return 1;
  585. }
  586.  
  587.  
  588. /* nextcwedge: return next clockwise edge from given edge around given face */
  589.  
  590. int nextcwedge (edge, face)
  591. int edge, face;
  592. {
  593.     switch (edge) {
  594.     case LB: return (face == L)? LF : BN;
  595.     case LT: return (face == L)? LN : TF;
  596.     case LN: return (face == L)? LB : TN;
  597.     case LF: return (face == L)? LT : BF;
  598.     case RB: return (face == R)? RN : BF;
  599.     case RT: return (face == R)? RF : TN;
  600.     case RN: return (face == R)? RT : BN;
  601.     case RF: return (face == R)? RB : TF;
  602.     case BN: return (face == B)? RB : LN;
  603.     case BF: return (face == B)? LB : RF;
  604.     case TN: return (face == T)? LT : RN;
  605.     case TF: return (face == T)? RT : LF;
  606.     }
  607. }
  608.  
  609.  
  610. /* otherface: return face adjoining edge that is not the given face */
  611.  
  612. int otherface (edge, face)
  613. int edge, face;
  614. {
  615.     int other = leftface[edge];
  616.     return face == other? rightface[edge] : other;
  617. }
  618.  
  619.  
  620. /* makecubetable: create the 256 entry table for cubical polygonization */
  621.  
  622. makecubetable ()
  623. {
  624.     int i, e, c, done[12], pos[8];
  625.     for (i = 0; i < 256; i++) {
  626.     for (e = 0; e < 12; e++) done[e] = 0;
  627.     for (c = 0; c < 8; c++) pos[c] = BIT(i, c);
  628.     for (e = 0; e < 12; e++)
  629.         if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
  630.         INTLIST *ints = 0;
  631.         INTLISTS *lists = (INTLISTS *) mycalloc(1, sizeof(INTLISTS));
  632.         int start = e, edge = e;
  633.         /* get face that is to right of edge from pos to neg corner: */
  634.         int face = pos[corner1[e]]? rightface[e] : leftface[e];
  635.         while (1) {
  636.             edge = nextcwedge(edge, face);
  637.             done[edge] = 1;
  638.             if (pos[corner1[edge]] != pos[corner2[edge]]) {
  639.             INTLIST *tmp = ints;
  640.             ints = (INTLIST *) mycalloc(1, sizeof(INTLIST));
  641.             ints->i = edge;
  642.             ints->next = tmp; /* add edge to head of list */
  643.             if (edge == start) break;
  644.             face = otherface(edge, face);
  645.             }
  646.         }
  647.         lists->list = ints; /* add ints to head of table entry */
  648.         lists->next = cubetable[i];
  649.         cubetable[i] = lists;
  650.         }
  651.     }
  652. }
  653.  
  654.  
  655. /**** Storage ****/
  656.  
  657.  
  658. /* mycalloc: return successful calloc or exit program */
  659.  
  660. char *mycalloc (nitems, nbytes)
  661. int nitems, nbytes;
  662. {
  663.    char *ptr = calloc(nitems, nbytes);
  664.    if (ptr != NULL) return ptr;
  665.    fprintf(stderr, "can't calloc %d bytes\n", nitems*nbytes);
  666.    exit(1);
  667. }
  668.  
  669.  
  670. /* setcenter: set (i,j,k) entry of table[]
  671.  * return 1 if already set; otherwise, set and return 0 */
  672.  
  673. int setcenter(table, i, j, k)
  674. CENTERLIST *table[];
  675. int i, j, k;
  676. {
  677.     int index = HASH(i, j, k);
  678.     CENTERLIST *new, *l, *q = table[index];
  679.     for (l = q; l != NULL; l = l->next)
  680.     if (l->i == i && l->j == j && l->k == k) return 1;
  681.     new = (CENTERLIST *) mycalloc(1, sizeof(CENTERLIST));
  682.     new->i = i; new->j = j; new->k = k; new->next = q;
  683.     table[index] = new;
  684.     return 0;
  685. }
  686.  
  687.  
  688. /* setedge: set vertex id for edge */
  689.  
  690. setedge (table, i1, j1, k1, i2, j2, k2, vid)
  691. EDGELIST *table[];
  692. int i1, j1, k1, i2, j2, k2, vid;
  693. {
  694.     unsigned int index;
  695.     EDGELIST *new;
  696.     if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
  697.     int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t;
  698.     }
  699.     index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
  700.     new = (EDGELIST *) mycalloc(1, sizeof(EDGELIST));
  701.     new->i1 = i1; new->j1 = j1; new->k1 = k1;
  702.     new->i2 = i2; new->j2 = j2; new->k2 = k2;
  703.     new->vid = vid;
  704.     new->next = table[index];
  705.     table[index] = new;
  706. }
  707.  
  708.  
  709. /* getedge: return vertex id for edge; return -1 if not set */
  710.  
  711. int getedge (table, i1, j1, k1, i2, j2, k2)
  712. EDGELIST *table[];
  713. int i1, j1, k1, i2, j2, k2;
  714. {
  715.     EDGELIST *q;
  716.     if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
  717.     int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t;
  718.     };
  719.     q = table[HASH(i1, j1, k1)+HASH(i2, j2, k2)];
  720.     for (; q != NULL; q = q->next)
  721.     if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
  722.         q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
  723.         return q->vid;
  724.     return -1;
  725. }
  726.  
  727.  
  728. /**** Vertices ****/
  729.  
  730.  
  731. /* vertid: return index for vertex on edge:
  732.  * c1->value and c2->value are presumed of different sign
  733.  * return saved index if any; else compute vertex and save */
  734.  
  735. int vertid (c1, c2, p)
  736. CORNER *c1, *c2;
  737. PROCESS *p;
  738. {
  739.     VERTEX v;
  740.     POINT a, b;
  741.     int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
  742.     if (vid != -1) return vid;                 /* previously computed */
  743.     a.x = c1->x; a.y = c1->y; a.z = c1->z;
  744.     b.x = c2->x; b.y = c2->y; b.z = c2->z;
  745.     converge(&a, &b, c1->value, p->function, &v.position); /* position */
  746.     vnormal(&v.position, p, &v.normal);               /* normal */
  747.     addtovertices(&p->vertices, v);               /* save vertex */
  748.     vid = p->vertices.count-1;
  749.     setedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
  750.     return vid;
  751. }
  752.  
  753.  
  754. /* addtovertices: add v to sequence of vertices */
  755.  
  756. addtovertices (vertices, v)
  757. VERTICES *vertices;
  758. VERTEX v;
  759. {
  760.     if (vertices->count == vertices->max) {
  761.     int i;
  762.     VERTEX *new;
  763.     vertices->max = vertices->count == 0 ? 10 : 2*vertices->count;
  764.     new = (VERTEX *) mycalloc((unsigned) vertices->max, sizeof(VERTEX));
  765.     for (i = 0; i < vertices->count; i++) new[i] = vertices->ptr[i];
  766.     if (vertices->ptr != NULL) free((char *) vertices->ptr);
  767.     vertices->ptr = new;
  768.     }
  769.     vertices->ptr[vertices->count++] = v;
  770. }
  771.  
  772.  
  773. /* vnormal: compute unit length surface normal at point */
  774.  
  775. vnormal (point, p, v)
  776. POINT *point, *v;
  777. PROCESS *p;
  778. {
  779.     double f = p->function(point->x, point->y, point->z);
  780.     v->x = p->function(point->x+p->delta, point->y, point->z)-f;
  781.     v->y = p->function(point->x, point->y+p->delta, point->z)-f;
  782.     v->z = p->function(point->x, point->y, point->z+p->delta)-f;
  783.     f = sqrt(v->x*v->x + v->y*v->y + v->z*v->z);
  784.     if (f != 0.0) {v->x /= f; v->y /= f; v->z /= f;}
  785. }
  786.  
  787.  
  788. /* converge: from two points of differing sign, converge to zero crossing */
  789.  
  790. converge (p1, p2, v, function, p)
  791. double v;
  792. double (*function)();
  793. POINT *p1, *p2, *p;
  794. {
  795.     int i = 0;
  796.     POINT pos, neg;
  797.     if (v < 0) {
  798.     pos.x = p2->x; pos.y = p2->y; pos.z = p2->z;
  799.     neg.x = p1->x; neg.y = p1->y; neg.z = p1->z;
  800.     }
  801.     else {
  802.     pos.x = p1->x; pos.y = p1->y; pos.z = p1->z;
  803.     neg.x = p2->x; neg.y = p2->y; neg.z = p2->z;
  804.     }
  805.     while (1) {
  806.     p->x = 0.5*(pos.x + neg.x);
  807.     p->y = 0.5*(pos.y + neg.y);
  808.     p->z = 0.5*(pos.z + neg.z);
  809.     if (i++ == RES) return;
  810.     if ((function(p->x, p->y, p->z)) > 0.0)
  811.          {pos.x = p->x; pos.y = p->y; pos.z = p->z;}
  812.     else {neg.x = p->x; neg.y = p->y; neg.z = p->z;}
  813.     }
  814. }
  815.